import os
import warnings
warnings.simplefilter("ignore")
import statsmodels.api as sm
from scipy import stats
import numpy as np
import pandas as pd
#
import matplotlib
import matplotlib.pyplot as plt
#
import cartopy.crs as ccrs
# 描画設定
plt.rcParams['figure.dpi'] = 300
plt.rcParams.update({'font.family': 'sans-serif', 'text.usetex': False,'pcolor.shading':'auto'})
def Mon():
import pandas as pd
mons = pd.DataFrame({
"name": ["Jan", "Feb", "Mar","Apr","May","Jun",\
"Jul","Aug","Sep","Oct","Nov","Dec","Jan"],
"days": [31,28,31,30,31,30,31,31,30,31,30,31,31],
"dayacc":[0,31,59,90,120,151,181,212,243,273,304,334,365],
"days_lp": [31,29,31,30,31,30,31,31,30,31,30,31,31],
"dayacc_lp":[0,31,60,91,121,152,182,213,244,274,305,335,366]
})
return mons
mons=Mon()
xticks=np.zeros(12)
xticks[:]=mons["dayacc"][:-1]
xlabels=mons["name"][:-1]
class Site:
def DomeF(self):
sx=14 ; sy=60 # 77°18′59″S 39°42′04″E
df_lat=-77.3
df_lon=39.66
return df_lon,df_lat,sx,sy
def DomeC(self):
sx=45 ; sy=59
dc_lat=-73.95
dc_lon=123.75
return dc_lon,dc_lat,sx,sy
dc_lon,dc_lat,dc_x,dc_y = Site().DomeC()
loadfile = "T2.npz"
dataset = np.load(loadfile)
T2 = dataset["T2"]
lon2 = dataset["lon2"]
lat2 = dataset["lat2"]
y = dataset["y"]
m = dataset["m"]
d = dataset["d"]
loadfile = "prcp.npz"
dataset = np.load(loadfile)
prcp = dataset["prcp"]
loadfile = "prcp_d18O.npz"
dataset = np.load(loadfile)
prcp_d18O = dataset["prcp_d18O"]
numofdays=366+365+365
model_dayst = len(y[y<2008])
model_dayed = model_dayst + numofdays
model_y = y [ model_dayst:model_dayed]
model_m = m [ model_dayst:model_dayed]
model_d = d [ model_dayst:model_dayed]
model_temp_df = T2 [dc_x-1, dc_y-1,model_dayst:model_dayed]
model_prcp_df = prcp [dc_x-1, dc_y-1,model_dayst:model_dayed]
model_prcp_d18O_df = prcp_d18O [dc_x-1, dc_y-1,model_dayst:model_dayed]
obs_DomeC = np.genfromtxt("obs_DomeC.csv",
delimiter=",", # 区切り文字
usecols=(0,1,2,3,4,5,6,7,8) # 読み込みたい列
)
obs_len = np.shape(obs_DomeC)[0]
obs_y_df = obs_DomeC[1:,0]
obs_m_df = obs_DomeC[1:,1]
obs_d_df = obs_DomeC[1:,2]
obs_temp_df = obs_DomeC[1:,3]
obs_prcp_df = obs_DomeC[1:,8]
obs_prcp_d18O_df = obs_DomeC[1:,5]
model_temp_dc2008=model_temp_df[model_y==2008]
model_prcp_dc2008=model_prcp_df[model_y==2008]
model_prcp_d18O_dc2008=model_prcp_d18O_df[model_y==2008]
obs_temp_dc2008=obs_temp_df[obs_y_df==2008]
obs_prcp_dc2008=obs_prcp_df[obs_y_df==2008]
obs_prcp_d18O_dc2008=obs_prcp_d18O_df[obs_y_df==2008]
model_temp_dc2009=model_temp_df[model_y==2009]
model_prcp_dc2009=model_prcp_df[model_y==2009]
model_prcp_d18O_dc2009=model_prcp_d18O_df[model_y==2009]
obs_temp_dc2009=obs_temp_df[obs_y_df==2009]
obs_prcp_dc2009=obs_prcp_df[obs_y_df==2009]
obs_prcp_d18O_dc2009=obs_prcp_d18O_df[obs_y_df==2009]
model_temp_dc2010=model_temp_df[model_y==2010]
model_prcp_dc2010=model_prcp_df[model_y==2010]
model_prcp_d18O_dc2010=model_prcp_d18O_df[model_y==2010]
obs_temp_dc2010=obs_temp_df[obs_y_df==2010]
obs_prcp_dc2010=obs_prcp_df[obs_y_df==2010]
obs_prcp_d18O_dc2010=obs_prcp_d18O_df[obs_y_df==2010]
def vert_lines(ax, xint):
if len(xint) < 2 :
ax.axvline(x=xint,ls='-',linewidth=1, color='k')
else:
for i in range(len(xint)):
ax.axvline(x=xint[i],ls='--',linewidth=1, color='k')
return ax
def hor_lines(ax, yint):
if len(yint) < 2 :
ax.axhline(y=yint,ls='-',linewidth=1, color='k')
else:
for i in range(len(yint)):
ax.axhline(y=yint[i],ls='--',linewidth=1, color='k')
return ax
fig = plt.figure(figsize=(15,10))
days=366
daylist=range(days)
col_obs = "black"
col_model = "red"
#-- d18Op --
ax = fig.add_subplot(3,1,1)
plt.title("(a)",loc="left",fontsize=20, weight="demibold")
## Model
ax.plot(daylist, model_prcp_d18O_dc2008, ms=10, color=col_model, linewidth=2, label="model")
ax.tick_params(labelsize=24)
ax.set_ylabel("$\mathsf{\delta^{18}O_p}$ modeled [\u2030]", fontsize=16, color=col_obs)
ax.set_xlim(0, days)
ax.set_ylim(-130,7)
ax.set_yticks([-120,-90,-60,-30,0])
ax.set_yticklabels([-120,-90,-60,-30,0],fontsize=16)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels,fontsize=16)
## Observation
ax2=ax.twinx()
ax2.scatter(daylist, obs_prcp_d18O_dc2008, color=col_obs, linewidth=2, alpha=0.5, label="observation")
ax2.tick_params(labelsize=24)
ax2.set_ylabel("$\mathsf{\delta^{18}O_p}$ observed [\u2030]", fontsize=16, color=col_obs)
ax2.set_xlim(0,days)
ax2.set_ylim(-80,-30)
ax2.set_yticks([-80,-60,-40])
ax2.set_yticklabels([-80,-60,-40],fontsize=16)
vert_lines(ax,xticks)
#-- Temprature --
ax = fig.add_subplot(3,1,2)
plt.title("(b)",loc="left",fontsize=20, weight="demibold")
## Model
ax.plot(daylist, model_temp_dc2008, ms=10, color=col_model, linewidth=2, label="model")
ax.tick_params(labelsize=24)
ax.set_ylabel("SAT modeled [\u00b0C]", fontsize=16, color=col_obs)
ax.set_xlim(0,days)
ax.set_ylim(-80,-10)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels,fontsize=16)
ax.set_yticks([-80,-60,-40,-20])
ax.set_yticklabels([-80,-60,-40,-20],fontsize=16)
## Observation
ax2=ax.twinx()
ax2.scatter(daylist, obs_temp_dc2008, color=col_obs, linewidth=2, alpha=0.5,label="observation")
ax2.tick_params(labelsize=24)
ax2.set_ylabel("SAT observed [\u00b0C]", fontsize=16, color=col_obs)
ax2.set_xlim(0, days)
ax2.set_ylim(-80,-10)
ax2.set_yticks([-80,-60,-40,-20])
ax2.set_yticklabels([-80,-60,-40,-20],fontsize=16)
vert_lines(ax,xticks)
#-- Precipitation --
ax = fig.add_subplot(3,1,3)
plt.title("(c)",loc="left",fontsize=20, weight="demibold")
## Model
ax.bar(daylist, model_prcp_dc2008, color=col_model, linewidth=2, alpha=0.8, label="model")
ax.tick_params(labelsize=24)
ax.set_ylabel("Precipitation\n[mm/day w.e.]", fontsize=16, color=col_obs)
ax.set_xlim(0,days)
ax.set_ylim(0,2.2)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels,fontsize=16)
ax.set_yticks([0,0.5,1,1.5,2])
ax.set_yticklabels([0,0.5,1,1.5,2],fontsize=16)
## Observation
ax2=ax.twinx()
ax2.bar(daylist, obs_prcp_dc2008, color=col_obs, linewidth=2, alpha=0.5, label="observation")
ax2.tick_params(labelsize=24)
#ax2.set_ylabel("Precipitation [mm/day w.e.]", fontsize=16, color=col_obs)
ax2.set_xlim(0,days)
ax2.set_ylim(0,2.2)
ax2.set_yticks([0,0.5,1,1.5,2])
ax2.set_yticklabels([0,0.5,1,1.5,2],fontsize=16)
vert_lines(ax,xticks)
ax.set_xlabel("2008", loc="center",fontsize=20)
h1, l1 = ax.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
plt.legend(h1+h2, l1+l2, loc='upper left',fontsize=20)
plt.subplots_adjust(top=0.92,bottom=0.1,left=0.1,right=0.9,hspace=0.3,wspace=0.05)
plt.show()
fig.savefig("Ant_proxy_DC2008.png")
fig = plt.figure(figsize=(15,10))
days=365
daylist=range(days)
col_obs = "black"
col_model = "red"
#-- d18Op --
ax = fig.add_subplot(3,1,1)
plt.title("(a)",loc="left",fontsize=20, weight="demibold")
## Model
ax.plot(daylist, model_prcp_d18O_dc2009, ms=10, color=col_model, linewidth=2, label="model")
ax.tick_params(labelsize=24)
ax.set_ylabel("$\mathsf{\delta^{18}O_p}$ modeled [\u2030]", fontsize=16, color=col_obs)
ax.set_xlim(0, days)
ax.set_ylim(-130,7)
ax.set_yticks([-120,-90,-60,-30,0])
ax.set_yticklabels([-120,-90,-60,-30,0],fontsize=16)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels,fontsize=16)
## Observation
ax2=ax.twinx()
ax2.scatter(daylist, obs_prcp_d18O_dc2009, color=col_obs, linewidth=2, alpha=0.5, label="observation")
ax2.tick_params(labelsize=24)
ax2.set_ylabel("$\mathsf{\delta^{18}O_p}$ observed [\u2030]", fontsize=16, color=col_obs)
ax2.set_xlim(0,days)
ax2.set_ylim(-80,-30)
ax2.set_yticks([-80,-60,-40])
ax2.set_yticklabels([-80,-60,-40],fontsize=16)
vert_lines(ax,xticks)
#-- Temprature --
ax = fig.add_subplot(3,1,2)
plt.title("(b)",loc="left",fontsize=20, weight="demibold")
## Model
ax.plot(daylist, model_temp_dc2009, ms=10, color=col_model, linewidth=2, label="model")
ax.tick_params(labelsize=24)
ax.set_ylabel("SAT modeled [\u00b0C]", fontsize=16, color=col_obs)
ax.set_xlim(0,days)
ax.set_ylim(-80,-10)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels,fontsize=16)
ax.set_yticks([-80,-60,-40,-20])
ax.set_yticklabels([-80,-60,-40,-20],fontsize=16)
## Observation
ax2=ax.twinx()
ax2.scatter(daylist, obs_temp_dc2009, color=col_obs, linewidth=2, alpha=0.5,label="observation")
ax2.tick_params(labelsize=24)
ax2.set_ylabel("SAT observed [\u00b0C]", fontsize=16, color=col_obs)
ax2.set_xlim(0, days)
ax2.set_ylim(-80,-10)
ax2.set_yticks([-80,-60,-40,-20])
ax2.set_yticklabels([-80,-60,-40,-20],fontsize=16)
vert_lines(ax,xticks)
#-- Precipitation --
ax = fig.add_subplot(3,1,3)
plt.title("(c)",loc="left",fontsize=20, weight="demibold")
## Model
ax.bar(daylist, model_prcp_dc2009, color=col_model, linewidth=2, alpha=0.8, label="model")
ax.tick_params(labelsize=24)
ax.set_ylabel("Precipitation\n[mm/day w.e.]", fontsize=16, color=col_obs)
ax.set_xlim(0,days)
ax.set_ylim(0,2.2)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels,fontsize=16)
ax.set_yticks([0,0.5,1,1.5,2])
ax.set_yticklabels([0,0.5,1,1.5,2],fontsize=16)
## Observation
ax2=ax.twinx()
ax2.bar(daylist, obs_prcp_dc2009, color=col_obs, linewidth=2, alpha=0.5, label="observation")
ax2.tick_params(labelsize=24)
#ax2.set_ylabel("Precipitation [mm/day w.e.]", fontsize=16, color=col_obs)
ax2.set_xlim(0,days)
ax2.set_ylim(0,2.2)
ax2.set_yticks([0,0.5,1,1.5,2])
ax2.set_yticklabels([0,0.5,1,1.5,2],fontsize=16)
vert_lines(ax,xticks)
ax.set_xlabel("2009", loc="center",fontsize=20)
h1, l1 = ax.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
plt.legend(h1+h2, l1+l2, loc='upper left',fontsize=20)
plt.subplots_adjust(top=0.92,bottom=0.1,left=0.1,right=0.9,hspace=0.3,wspace=0.05)
plt.show()
fig.savefig("Ant_proxy_DC2009.png")
fig = plt.figure(figsize=(15,10))
days=365
daylist=range(days)
col_obs = "black"
col_model = "red"
#-- d18Op --
ax = fig.add_subplot(3,1,1)
plt.title("(a)",loc="left",fontsize=20, weight="demibold")
## Model
ax.plot(daylist, model_prcp_d18O_dc2010, ms=10, color=col_model, linewidth=2, label="model")
ax.tick_params(labelsize=24)
ax.set_ylabel("$\mathsf{\delta^{18}O_p}$ modeled [\u2030]", fontsize=16, color=col_obs)
ax.set_xlim(0, days)
ax.set_ylim(-130,7)
ax.set_yticks([-120,-90,-60,-30,0])
ax.set_yticklabels([-120,-90,-60,-30,0],fontsize=16)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels,fontsize=16)
## Observation
ax2=ax.twinx()
ax2.scatter(daylist, obs_prcp_d18O_dc2010, color=col_obs, linewidth=2, alpha=0.5, label="observation")
ax2.tick_params(labelsize=24)
ax2.set_ylabel("$\mathsf{\delta^{18}O_p}$ observed [\u2030]", fontsize=16, color=col_obs)
ax2.set_xlim(0,days)
ax2.set_ylim(-80,-30)
ax2.set_yticks([-80,-60,-40])
ax2.set_yticklabels([-80,-60,-40],fontsize=16)
vert_lines(ax,xticks)
#-- Temprature --
ax = fig.add_subplot(3,1,2)
plt.title("(b)",loc="left",fontsize=20, weight="demibold")
## Model
ax.plot(daylist, model_temp_dc2010, ms=10, color=col_model, linewidth=2, label="model")
ax.tick_params(labelsize=24)
ax.set_ylabel("SAT modeled [\u00b0C]", fontsize=16, color=col_obs)
ax.set_xlim(0,days)
ax.set_ylim(-80,-10)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels,fontsize=16)
ax.set_yticks([-80,-60,-40,-20])
ax.set_yticklabels([-80,-60,-40,-20],fontsize=16)
## Observation
ax2=ax.twinx()
ax2.scatter(daylist, obs_temp_dc2010, color=col_obs, linewidth=2, alpha=0.5,label="observation")
ax2.tick_params(labelsize=24)
ax2.set_ylabel("SAT observed [\u00b0C]", fontsize=16, color=col_obs)
ax2.set_xlim(0, days)
ax2.set_ylim(-80,-10)
ax2.set_yticks([-80,-60,-40,-20])
ax2.set_yticklabels([-80,-60,-40,-20],fontsize=16)
vert_lines(ax,xticks)
#-- Precipitation --
ax = fig.add_subplot(3,1,3)
plt.title("(c)",loc="left",fontsize=20, weight="demibold")
## Model
ax.bar(daylist, model_prcp_dc2010, color=col_model, linewidth=2, alpha=0.8, label="model")
ax.tick_params(labelsize=24)
ax.set_ylabel("Precipitation\n[mm/day w.e.]", fontsize=16, color=col_obs)
ax.set_xlim(0,days)
ax.set_ylim(0,2.2)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels,fontsize=16)
ax.set_yticks([0,0.5,1,1.5,2])
ax.set_yticklabels([0,0.5,1,1.5,2],fontsize=16)
## Observation
ax2=ax.twinx()
ax2.bar(daylist, obs_prcp_dc2010, color=col_obs, linewidth=2, alpha=0.5, label="observation")
ax2.tick_params(labelsize=24)
#ax2.set_ylabel("Precipitation [mm/day w.e.]", fontsize=16, color=col_obs)
ax2.set_xlim(0,days)
ax2.set_ylim(0,2.2)
ax2.set_yticks([0,0.5,1,1.5,2])
ax2.set_yticklabels([0,0.5,1,1.5,2],fontsize=16)
vert_lines(ax,xticks)
ax.set_xlabel("2010", loc="center",fontsize=20)
h1, l1 = ax.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
plt.legend(h1+h2, l1+l2, loc='upper left',fontsize=20)
plt.subplots_adjust(top=0.92,bottom=0.1,left=0.1,right=0.9,hspace=0.3,wspace=0.05)
plt.show()
fig.savefig("Ant_proxy_DC2010.png")
import statsmodels.tools.eval_measures as smte
def corr(a1,a2):
aa1=a1[~(np.isnan(a1)+np.isnan(a2))]
aa2=a2[~(np.isnan(a1)+np.isnan(a2))]
return stats.spearmanr(aa1,aa2)[0]
obs_prcp_d18O_dc=np.concatenate([obs_prcp_d18O_dc2008,obs_prcp_d18O_dc2009,obs_prcp_d18O_dc2010])
model_prcp_d18O_dc= np.concatenate([model_prcp_d18O_dc2008,model_prcp_d18O_dc2009,model_prcp_d18O_dc2010])
obs_prcp_dc=np.concatenate([obs_prcp_dc2008,obs_prcp_dc2009,obs_prcp_dc2010])
model_prcp_dc= np.concatenate([model_prcp_dc2008,model_prcp_dc2009,model_prcp_dc2010])
obs_temp_dc=np.concatenate([obs_temp_dc2008,obs_temp_dc2009,obs_temp_dc2010])
model_temp_dc= np.concatenate([model_temp_dc2008,model_temp_dc2009,model_temp_dc2010])
corr_d18O = "{:4.2f}".format(corr(obs_prcp_d18O_dc,model_prcp_d18O_dc))
corr_temp = "{:4.2f}".format(corr(obs_temp_dc,model_temp_dc))
corr_prcp = "{:4.2f}".format(corr(obs_prcp_dc,model_prcp_dc))
df=pd.DataFrame.from_dict({
"$r_{\mathsf{\delta^{18}O_p}}$":corr_d18O,
"$r_{\mathsf{SAT}}$" :corr_temp,
"$r_{\mathsf{prcp}}$" :corr_prcp,
},orient="index")
df.columns=["Correlation coefficients (model vs. observation)"]
df
def sqrt(a1,a2):
aa1=a1[~(np.isnan(a1)+np.isnan(a2))]
aa2=a2[~(np.isnan(a1)+np.isnan(a2))]
return smte.rmse(aa1,aa2)
df=pd.DataFrame.from_dict({
"$\mathsf{\delta^{18}O_p}$":[sqrt(obs_prcp_d18O_dc,model_prcp_d18O_dc),
np.nanmean(obs_prcp_d18O_dc),np.nanmean(model_prcp_d18O_dc),
np.nanmedian(obs_prcp_d18O_dc),np.nanmedian(model_prcp_d18O_dc),
np.nanstd(obs_prcp_d18O_dc),np.nanstd(model_prcp_d18O_dc)],
"$\mathsf{SAT}$" :[sqrt(obs_temp_dc,model_temp_dc),
np.nanmean(obs_temp_dc),np.nanmean(model_temp_dc),
np.nanmedian(obs_temp_dc),np.nanmedian(model_temp_dc),
np.nanstd(obs_temp_dc),np.nanstd(model_temp_dc)],
"$\mathsf{prcp}$" :[sqrt(obs_prcp_dc,model_prcp_dc),
np.nanmean(obs_prcp_dc),np.mean(model_prcp_dc),
np.nanmedian(obs_prcp_dc),np.median(model_prcp_dc),
np.nanstd(obs_prcp_dc),np.nanstd(model_prcp_df)],
},orient="index")
df.columns=["Sqrt", "mean(obs)","mean(model)","median(obs)","median(model)","SD (obs)","SD(model)"]
df
Kanon Kino (kanon@aori.u-tokyo.ac.jp)